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A METHOD FOR OBTAINING ANALYTICAL SOLUTIONS TO THE EQUATION 
FOR WIND-DRIVEN CIRCULATION IN A SHALLOW SEA OR LAKE 
by Marvin E. Goldstein, Willis H. Braun, and Richard T. Gedney 

Lewis Research Center 

SUMMARY 

A technique for obtaining solutions to a certain class of partial differential equations 
is introduced. This class of equations includes Welander T s equation for the wind-driven 
circulation in shallow seas and lakes for a large class of bottom topographies. The 
technique, together with some results based on conformal mapping, is used to reduce 
the problems of finding solutions to Welander T s equation for a closed body of water to the 
problem of solving an ordinary differential equation. To illustrate the principles in- 
volved, the method is used to obtain the complete analytical solution for an elliptically 
shaped body of water with a particular choice of bottom topography when the depth is 
greater than one-half the Ekman friction depth. 


INTRODUCTION 

Lake currents are very important in determining what chemical and biological 
changes will take place in a lake, since they control the distributions of sediment and pol- 
lution which have been introduced into a lake. However, it is a major and expensive task 
to measure the lake currents in large lakes at very many locations. It is, therefore, 
important to be able to predict these currents analytically. Welander (ref. 1) derived 
the basic partial differential equation which describes the wind-driven circulation in a 
shallow sea or lake. 

In this report a method is developed for obtaining product solutions to a certain 
class of partial differential equations. It is shown that, if the function which describes 
the distribution of depth of a lake or sea belongs to a certain quite general class of 
functions, then the technique developed herein can be used to reduce the problem of 
solving the partial differential equation and boundary conditions which describe the wind- 
driven circulation in the lake or sea to the problem of solving an ordinary differential 



equation and its boundary conditions. To illustrate the principles involved, the technique 
is used to obtain the complete analytical solution for an elliptically shaped body of water 
with a particular choice of bottom topography when the depth is greater than half the 
Ekman friction depth. 

The report begins by considering the general canonical form of the second-order 
linear elliptic partial differential equation in two independent variables. It is shown that 
if the two combinations of coefficients of the equation which are known as the invariants 
have a certain functional form then the equation can be transformed into an equation 
whose coefficients are functions of only one of its independent variables. Equations of 
this type, although they are not necessarily separable, possess infinitely many solutions 
which are products of a function of one of the variables with a function of the other. 

These solutions can be found by solving a single ordinary differential equation. Thus, an 
infinite family of solutions to the original differential equation can be found. It is antici- 
pated that this family is sufficiently rich to be complete. However, instead of pursuing 
these matters further in this report, a new transformation of the independent variables 
is introduced. The transformed equation is again one whose coefficients depend on only 
one of the independent variables. A quite general boundary value problem for this equa- 
tion is then posed in the interior of the unit circle in the plane for which the new variables 
are polar coordinates. K is then shown that by taking finite Fourier transforms the 
problem can be reduced to the problem of solving a second-order ordinary differential 
equation. 

Next, the partial differential equation which describes the wind-driven circulation in 
a shallow sea or lake is introduced. A simplified form of this equation which is valid 
when the depth of the body of water is larger than half the Ekman friction depth is also 
given. It is shown that whenever the bottom topography of the body of water can be de- 
scribed by a function which belongs to a certain general class of functions these equations 
are of the type discussed herein. It is also shown that if, in addition, the depth of the 
lake is constant along the shore then the boundary value problem for determining the 
wind-driven circulation in the body of water can be transformed into a boundary value 
problem of the type discussed above and can, therefore, be reduced to the problem of 
solving an ordinary differential equation. The coefficients in this equation depend on the 
bottom topography. The ordinary differential equation together with its boundary condi- 
tions are solved for a particular choice of the bottom topography. Finally, the complete 
solution to the problem is obtained for an elliptically shaped lake or sea. The transfor- 
mations which are carried out to obtain the solutions are best accomplished by using 
certain techniques of conformal mapping. 



SYMBOLS 


function of h/d given by eq. (42) 

function of x, y 

defined in eq. (4) 

function of h/d given by eq. (42) 

function of x, y 

defined in eq. (4) 

function of h/d defined in eq. (37) 
boundary of sea or lake 
function of x, y 
defined in eq. (4) 

function of h/d defined in eq. (37) 

Ekman friction depth, ir(2v/t c )^^ 
functions of h/d defined in eq. (38) 
defined by eq. (17) 
function of x, y 
Coriolis parameter 
acceleration due to gravity 
defined in eq. (49) 

value of H at lake or sea boundary (a constant) 
depth of lake or sea 
defined in eq. (114) 

invariant of linear, elliptic partial differential equation, defined in eq. 
Bessel function of first kind of order n 

invariant of linear, elliptic partial differential equation, defined in eq. 

complete elliptic integral of first kind 

K(k') 

kernel in solution for £2 n (r), eq. (104) 
kernel in solution for ^(r), eq. (94) 



k modulus of elliptic integral 


k* 
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p,q 

Qp ^2 

q P q 2 
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s 

t 

U 


u 

V 

v 

W 


w 


X 

Y 

y 

z 


z 


semimajor and semiminof axes of ellipse 
constants of integration 

constants, eq. (96) 

constant 

integer 

unit normal to sea or lake boundary 
x- and y- components of n 
constants 

x- and y-eomponents of volume flow 

x- and y- components of velocity 

radial distance in x~Pl ane 

X r 
n 

dummy variable 

length of curve in z-plane 

constant 

complex variable, § + ir] 
dependent variable 
harmonic function of x, y 
dependent variable 
harmonic conjugate to u 
Wronskian 

complex function, u + iv 

independent variable, coordinate of lake-surface plane 
Bessel function of second kind of order n 
independent variable, coordinate of lake-surface plane 
vertical coordinate 
x + iy 
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z 0 
a, p 

r n 
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? 
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n 




a 

(1) (2) 
ex' . cr ’ 


t 1’ t 2 


T t T t 

T l> r 2 
<P 
X 
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n 

oP 

“n 

u> 

J 2 ) 

n » n 


a particular point interior to 9d 
defined in eq. (51) 

Fourier coefficient of &(r , 0), eq. (78) 

periodic function of 6 

Fourier coefficient of y(6) 

parameter in depth distribution, eq. (84) 

displacement of lake or sea surface 

imaginary part of t 

angular coordinate in x _ P lane 

wind direction 

complex constant 

defined in eq. (95) 

defined in eq. (128) 

coefficient of vertical eddy viscosity 

real part of t 

radial polar coordinate 

„<« . i„< 2 > 

linear combinations of wind stress, eq. (36) 
x- and y- components of wind stress 

magnitude of dimensionless wind stress 

27TTj/gd, 27TT 2 /gd 
differentiable function of u 
defined in eq. (20) 
function of u 

Fourier coefficient of V, eq. (27) 

complementary function for £2 n 
particular solution of differential equation for £2^ 
function of x, y 
defined by eq. (102) 
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Subscripts: 

r.O 


u, v V 


partial derivatives 


x,y 


J 


Superscript: 


denotes complex conjugate 


SOLUTIONS OF A CERTAIN GENERAL CLASS OF 
PARTIAL DIFFERENTIAL EQUATIONS 


The most general, second order, linear, elliptic, partial differential equation in two 
independent variables can always be reduced to the canonical form 


i!» + 3?M +a lU +b 3U +cU = f 


dx 


dy 


dx 


dy 


( 1 ) 


where the coefficients a, b, and c and the nonhomogeneous term f are functions of x 
and y. 

Any change of the dependent variable of the form 

V = e w U (2) 

where oo is any sufficiently differentiable function of x and y transforms equation (1) 
into an equation of the same type. Thus, the function V satisfies the differential equa- 
tion 


Ar*aiY + a.W + b" + 5v«e"f 


dx‘ 


ay 


3x 


3y 


( 3 ) 


where the coefficients a, b, and c are related to the coefficients a, b, and c of the 
original equation (1) by 
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a = a - 2a> x 

b = b - 2w y > 

c = c - + Wyy) - (aco x + bo> y ) + (o> x ) 2 + (w y ) 2 

It is easy to show (ref. 2) by eliminating the function oj between these relations 

that 


(4) 


a - b ~ a - b 
y x y x 


C - - (2 + b ) - ±(2 2 + b 2 ) = c - i(a v + bj - — (a 2 + b 2 ) 


Thus, the quantities $ and $ defined in terms of the coefficients of any equation of 
the form equation (1) by 


* - a y - \ 


(5) 




( 6 ) 


are invariant under a transformation of the type shown in equation (2). They are, there- 
fore, called the invariants of the differential equation (1) (ref. 3, vol. VI, art. 191). 

We shall now restrict the discussion to the class of all equations of the type in equa- 
tion (1) whose invariants and £ can be expressed in the form 

& = <p’(u) [(u x ) 2 + (u y ) 2 ] (7) 


/ = *(U) [(u x ) 2 + (U y ) 2 ] (8) 

where u is any nonconstant harmonic function of x and y (i. e. , u^ + u = 0), \j/ is 
any function of u and <p is any differentiable function of u. Since u satisfies 
Laplace 1 s equation, it is easy to see from definition (5) and equation (7) that 


i i t a - vH -s;[ b + vw] 


9x 2 1 


This shows that there exists a differentiable function w of x and y which is deter - 
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mined to within an unimportant constant by 


< j ) 


x 


CO 


y 


j[ a - vH 

i[b * V <»)] 


'\ 




(9) 


Hence, when the change of variable (2), with w determined from equations (9), is ap- 
plied to equation (1) it follows that the coefficients a, b, and c of the resulting equa- 
tion (3) are given by 


a = u y <p(u) 


b = - u x^ u) 


c = 


[,(U) + ^<»(U)J [(U/ + (Uy)^ 


Hence, V satisfies the equation 


V + V + <p(u)(u V - u V ) + 
xx yy ^ v ; ' y x x y ; 


*(u) + 


i^(u,] f 


(U X ) 2 + (U y ) 2 


„ o>. 

V = e f 


( 10 ) 


Now let v be the harmonic conjugate* to u. Then u and v satisfy the Cauchy- 
Riemann equations 


u = v 
x y 


u = - V 

y x 


(ID 


It follows from these equations that the Jacobian 9(u, v)/8(x, y) is given by 


9(u, y) = 
3(x, y) 


u x u y 
v x v y 


= (U,) 2 + ^y) 2 


(12) 


*Thus v satisfies Laplace's equation and the function w = u + iv is an analytic 
function of the complex variable z = x + iy. 
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Since by hypothesis u is nonconstant, it follows from equation (12) that 


3(x,y) 


Hence, equation (10) can be further transformed by changing the independent vari- 
ables x and y to the independent variables u and v. Now it is well known (ref. 4, 
p. 187) that under this change of variable the Laplacean of V transforms as 2 

V xx + V yy = < V uu + V vv> [<“/ * <V *] < 13 > 

Since the chain rule implies that 


V = V u + V v 
x u x v x 


V = V u + V v 
y u y v y 


it follows from the Cauchy-Riemann equations (11) that 

u V - u V = -V f"(u ) 2 + (u ) 2 1 
y x x y vL' x' v y' J 


It is easy to show that^ 


(14) 


2 

The convention of writing V(u, v) in place of V(x(u, v), y(u, v)) will be adopted since 
no confusion is likely to result from this. 

3 

Equations (11) imply that x and y also satisfy the Cauchy-Riemann equations with 
respect to u and v. Thus, = y v and y u = -x y . Hence, 


But 


Hence, 


dfoy) = 

a(u, v) 


(V 2 + (*v )2 


d(u, v) 9(x, y) = ! 

3(X, y) a(u, v) 


(u x ) 2 - (V 2 


1 

(V 2 + (V 2 
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1 


(15) 


(V 2 + (x v )2 = 


( u x ) 2 + ( U y ) 2 


Upon substituting equations (13) to (15) into equation (10) we find, since 


2 2 

(u x ) + (u y ) ^ 0, that V must satisfy the equation 


V uu + V vv - ^ u > V v + I + 


|V(u) +i<p 2 (u) V = 


§ 


where 


#■= [<\) 2 + < V 2 l 


w. 
e f 


(16a) 


(16b) 


Notice that the coefficients of this equation are functions only of the independent 
variable u and not of the independent variable v. Then direct substitution shows that, 
for every real or complex constant X, the associated homogeneous equation (i. e. , 
eq. (16a) with & ■=■ 0) possesses a solution of the form 


V(u, v) = e Xv U x (u) 


(17) 


provided that satisfies the ordinary differential equation 

2 




(*-!«,) + *Ju x = o 


Of course, any linear combination of these "elementary" solutions is also a solu- 
tion. However, the question of whether these elementary solutions form a complete set 
of solutions (i. e. , whether every solution to eq. (16a) with ^=0 can be expressed as a 
linear combination of these elementary solutions) is beyond the scope of this report. The 
use of linear combinations of elementary solutions of the type shown in equation (17) to 
obtain solutions to specific boundary value problems involving homogeneous equations 
whose coefficients depend only on one of the independent variables is illustrated in ref- 
erences 5 and 6. Instead of pursuing these matters further, we shall turn to a discussion 
of an important special type of boundary value problem for the nonhomogeneous equation 
(16a). Before this is done, however, it will be convenient to introduce one further 
change of variable. 
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SOLUTIONS OF A CERTAIN TYPE OF BOUNDARY VALUE PROBLEM 


To this end notice that, since u and v are conjugate harmonic functions, the 
function 


w = u + iv (18) 

is an analytic function of the complex variable 

z = x + iy (19) 

Hence, the function x defined by 

X = e w (20) 

is also an analytic function of the complex variable z. Expressing x * n the polar form 

X=re i9 (21) 

reveals that 

u = In r 1 

> (22) 

v = 9 J 

When the change of variables defined by equation (22) is introduced into equation 
(16a), we obtain (where the convention introduced previously of writing cp(r) in place of 
cp(ln r) etc. is used) 

r — (r — ) + — - - cp( r) — + T^(r) + --(p 2 (r)l V = 3<{y, 9) (23) 

9r \ 3r / d9 [_ 4 J 

We shall now consider the following boundary value problem for equation (23) in the 
unit circle of the x~Pl aile shown in figure 1. 
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Figure 1. -x-plane. 


r ( rV r)r + V ee~ (p{r)v d + 


V(r) 

4 


i^(r) + -cp Z ( r) V = 


0 < r < 1 


- 77 < 9 < 77 


pV r + qVg + sV = y(0) for r = 1, -tt < 0 < 77 


(24) 


(25) 


where p, q, and s are constants and y is a periodic function of 6. We require that 
V be a continuously differentiable function within the unit circle. Thus, in particular 
V(0, 6) must be finite and independent of 0 , and the following periodicity conditions must 
hold: 


V(r, -77) = V(r, 77) 
V 0 (r, -77) = V 0 (r, 77) 


0 < r < 1 


(26) 


In order to obtain a solution to this problem put 
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o (r) = — r V(r, 0)e _in0 d0 n = 0, ±1, ±2, . . . (27) 

n 2 tt J-k 

Then it follows from the theory of Fourier series that this transform can be inverted to 
obtain 


V(r, Q) = ^ syr)e in0 
n=-°° 


(28) 


Upon integrating by parts, we find that 


— f V fl (r, 9)e~ ind dd = [v(r, n) - V(r, -7r)l + inft (r) 

2tt y 2 tt l j 


Hence, the first periodicity condition in equation (26) shows that 


-- H V 0 (r, 0) e _in0 d0 = inn n 

27T •'-7T 


(r) 


(29) 


In a similar way it follows from the second periodicity condition in equation (26) that 


r- /"v ra (r,9)e- tae de=-n 2 syr) 

*S7T % '-TT 


(30) 


We shall now show that each function O n can be determined as the solution of a 
certain ordinary differential equation. To this end, multiply equations (24) and (25) by 
e”in^/27r and integrate both sides of the resulting expressions with respect to 0 be- 
tween -it and t t. We find, upon making the definitions 


r n (r) s — r ® (r > e)e_in0 

n 2 t t J-n 


d0 


> n = 0, ±1, ±2, . . . 


- 1 f n /n\„-in0 
y = — / y(0)e 

11 27 T J-TT 


d0 


(31) 


and using equations (29) and (30), that must satisfy the ordinary differential equation 
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n = 0, ±1, ±2, . . . 


(32) 


r(r^)' + +~(p 2 - n 2 - in yj fi n = T n 

subject to the boundary condition 

pST (1) + (s + inq)S2 n (l) = r n n = 0, ±1, ±2, . . . (33) 

Equation (32) and the single boundary condition in equation (33) are not sufficient to 
completely determine Notice, however, that the origin (r = 0) is a regular singular 

point of the ordinary differential equation (32) (we are, of course, assuming that the 
functions cp and i p are bounded at r = 0). Hence, one of the homogeneous solutions to 
this equation will be bounded at r = 0 while the other will be unbounded (ref. 7, pp. 101- 
108). But since V must be bounded at r = 0, each £2 n must also be bounded. Hence, 
there must exist (for each integer n) a finite number M n such that 

lim £2 n (r) < M n n = 0, ±1, ±2, . . . (34) 

r— 0 


and the solution to the nonhomogeneous equation (32) which satisfies the boundedness 
condition in equation (34) will involve only one arbitrary constant. This constant is 
determined by the boundary condition shown in equation (33). Thus, the Fourier coeffi- 
cients S2 n are completely determined by the differential equations (32) together with the 
conditions in equations (33) and (34). When these solutions for n = 0, ±1, ±2, . . . 
have been found, they can be substituted into equation (28) to yield the complete solution 
to the boundary value problem posed by equations (24) and (25). 


BASIC EQUATIONS FOR WIND-DRIVEN CIRCULATION IN SHALLOW SEA OR LAKE 

The basic equation governing the wind-driven circulation in a one-layer shallow sea 
or lake was derived by Welander (ref. 1). The basic assumptions used in deriving the 
equation from the Navier-Stokes equations are that the water density is constant, the 
vertical eddy viscosity is independent of depth but dependent on wind velocity, the pres- 
sure is hydrostatic, and the lateral friction and nonlinear acceleration terms can be 
neglected. A detailed discussion of the derivation is contained in reference 8. Welander 
wrote the equation with the surface displacement as the dependent variable. This equa- 
tion for the surface displacement £ of a lake or sea of variable depth h is (ref. 1) 
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(ghE)(? xx + ^yy) + [fe hE >x + (S hF )y]? x + [fehE) y - (ghF) x ]? y 

= - fa M -faa< 2 >! (35) 

l2n J x |_2ir I 


where x and y are the ho rizontal, locally Cartesian coordinates, g is the acceleration 
due to gravity, d = ny2u /f c is the Ekman friction depth, v is the coefficient of vertical 
eddy viscosity, and f c is the Coriolis parameter which introduces the effect of the 
Earth's rotation. In addition, 


t (D, 


27TTj 

gd 


C - 


2lJT 2 

gd 




D 


y (2) = 2lTT l 

gd 


D + 


2ttt2 

gd 


(36) 


where Tj and 


T2 are, respectively, the x- and y-components of the wind stress and 



1 

sin ^27 r — ^ - sinh ^277 

1 *\ 

(277 
V dy 

| cos ^ + cosh ^2 tt — ^ 


1 

sin ^27r — ^ + sinh ^27r — ^ 

_i_ 1 

H): 

cos ^2 tt + cosh ^2tt — ^ 

T L 


(37) 


(38) 
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The x- and y-components of the total volume flow and Qg, respectively, are 
given in terms of the surface displacement in reference 1 as 


Q 1 H r q l dZ H — T a(1) + — (E? x " F vl 
1 *^-h L 2irf L d x y I 

Q 2 -/, 0 ^ dZ -^P 2)+ f (F? ^ E v] 


( 39 ) 


The boundary condition for equation (35) is obtained by requiring that the normal volume 
flow vanish at the coasts. Thus, for any closed lake or sea let n^ and ng denote the 
x and y components, respectively, of the outward drawn unit normal n to its boundary. 
Then the boundary condition for equation (35) is 


Q l n l + Q 2 n 2 = 0 


(40) 


at all points of the boundary of the body of water. 

It is convenient to rearrange equation (35) to obtain 


p + if + 

>xx Vy 


— 



" 


A (277 — ^ - B 



a l 2 , a) 

+ B ^2t7 — \ 

' d 4 

l 

X 

\ d). 

, v d 4 


? y = f 


(41) 


where 


2 sin 


A s 


27 r -^sinh^27T — ^ 


cosh (2tt — \ + cos (27 7 — 


sinh ^277 


sHH) 




B 5 


sinh ( 277 - \ + sin /2 tt h 


COSh (277 - \ + cos ^277 — 


f [4 11 + 4 2) ] 


277 — E 

d 


(42) 


(43) 
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It can be seen from this figure that for h/d >1/2 


2tt - F * 2tt - - 1 
d d 


2tt — E * -1 
d 

Thus, for h/d >1/2 expression (43) for f and equation (39) for the components of the 
total volume flow can be approximated fairly closely by the simpler expressions 

f » o (1) + a< 2 > for £ >± (45) 

x y d 2 


and 


for ~>± (46) 

d 2 


Q x ~ — - 
2 rf„ 


k (1) - ? v + fl - 2 tt% 


Q, 


a gd 

27rf 


t ( 2 ) 


t - 1-2^ -k 


The coefficients C and D of equation (37) are plotted in figure 4. It is also shown in 
this figure that for h/d >1/2 the asymptotic forms 

C - 2e" 7rh/d sin ^ 
d 

D « 2e _7rh//d cos — - 1 
d 


agree with the exact expressions for C and D to a fairly close approximation. Hence, 

ad) , 

2e~ 7Th/d 


equation (36) for ad) and a^ 2 ) can be approximated by 


T d) 27rT i 


sin — + 


2^ TT, 


2 (l - 2e -7rh/d cos — 


gd 

y (2) a 27fT l 

gd 


gd 


2e 


■irh/d CQS nh 


2tt T t 


/ for — > i 


1 + 


gd 


2 2 e- ,rh/d sin I* 


(47) 
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Figure 4. - Coefficients C and D. 


It can also be seen from figure 4 that, for h/d >1, C and D can be further approxi 
mated by the limits for h/d — namely, 

C *0 

and 

D ~ -1 

with a fair degree of accuracy. Hence, equation (36) can be replaced by the even sim 
pier expressions 


( 48 ) 


(1) 2 tit 2 ''l 

a * ^ 

gd 

>• 

(2) 2t rr x 

a « - 


gd J 



CLASS OF DEPTH DISTRIBUTIONS FOR WHICH EQUATIONS CAN BE SOLVED 

BY FINITE FOURIER TRANSFORMS 


We shall now suppose that there exists a nonconstant harmonic function u of x and 
y and an arbitrary function H of u only such that the depth distribution h can be ex- 
pressed in the form 


= _L H[u(x,y)] (49) 

d 2u 

This is a fairly general functional form and it will be possible, for any one of a large 
number of lakes and seas, to choose the functions u and H in equation (49) so that the 
depth distribution is approximated fairly closely by a relation of this type. 

When the relation (49) is substituted into equations (41) and (42), we find that 


^xx + ^yy + 2 [ a < u > u x " £( u ) u y ? x + 2 [«( u ) u y + 0 ( u ) u x]^y = f (50) 


where 


sin H sinh H dH 

(cosh H + cos H)(sinh H - sin H) du 

► 

0 _ _1 sinh H + sin H dH 
2 cosh H + cos H du 


j = cosh H + cos H 
sinh H - sin H 


9x 9y 


(51) 


(52) 


and it follows from the preceding discussion and equations (44) that, when h/d > 1/2, the 
coefficients a and /3 are given to a close approximation by 
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a ~ 0 


2 du 



( 53 ) 


Now the invariants of equation (50) satisfy conditions (7) and (8) for all differentiable 
functions a and /3. Hence, this equation is of the type discussed in the preceding sec- 
tion. In order to show this, notice that since 


9y 


ja (u) u x ] = ~ [a(u)u y ] 


it follows from definitions (5) and (6) that the invariants $ and £ of equation (50) are 




/ = [«W“x] + [“< u)u y]} - {H] 2 + + <V 2 ] 

Or, since u satisfies Laplace’s equation, 


- 2 S? K + <v 2 ] 

^ = -(5f +a2 ^)N 2 + < V 2 ] 


(54) 

(55) 


Thus, the invariants & and of equation (50) satisfy conditions (7) and (8) with 

cp(u) = -20(u) (56) 

^(u) = - {[*(*)]' + Jo!(u)J 2 + [/3( u )] 2 j (57) 

Equation (50) can, therefore, be transformed into an equation of the type (16a) and, 
hence, also into an equation of type (23). Thus, it follows from equations (1), (9), (50), 
and (56) that 


w 


X 


era 


x 
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CO = 

y 


QO 1 


Hence, 


CO = 



du 


Therefore, definition (51) shows that 




sin H sinh H 


(cosh H + cos H)(sinh H - sin H) 


dH 


or 



sinh H - sin H 
cosh H + cos H 



and equation (38) and figure 3 show that 


(58) 


u> ~ 0 if — > — 


(59) 


Therefore, in view of equations (1), (2), (38), (41), and (49) 


V = 


V sinh H - sin H 
cos H + cosh H 


and 


V » ? 


if *>I 


d 2 

The introduction of equations (56) and (57) into equation (17a) yields 


V uu + V vv - 2 ^< U)V V 


where, from equations (16b), (52), and (58), 


- |[a(u) ' + a(u)J 2 j>V = 


9 


9 - [(x u ) 2 + (x) 2 ] /c° shH+ c osH pi) + J2)1 

L u v j y gjjjjj jj + g h L x y J 


(60) 


(61) 


(62) 


(63) 


23 



iiii min ii ii 


mu inn lllll■■llll■lll mi 


i mu i i n ii 


Upon using the polar coordinates r and 0, equation (23) specializes to 

r 1 / r w] t * t Wr) iV f r Wr) t U)l4 = ji (64) 

dr \ dr / dQ 2 dd [ dr L J J 


When h/d > 1/2, it follows from equations (53), (61), and (22), and figure 3 that this 
equation becomes approximately 


r i/ r &) + A + r !? = / 

dr \ dr/ g0 2 dr dd 


9= 


(x u ) 2 + (x v ) 2 ][ 4 1) 




(65) 


We shall now consider a lake or shallow sea bounded by a closed curve such as il- 
lustrated in figure 5. We shall suppose that the bottom topography of the basin can be 
approximately described by a function of the form (49) and that the depth along the shore 
is constant. Hence, the harmonic function u is constant along r £ . Now let v be 


y 
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the harmonic conjugate to u and put 


w = u + iv 

Then w is an analytic function of z = x + iy, and the function x defined by equation (20) 
is also an analytic function of z. But 


u 


X = e 


Hence, 


| x(z) | = constant for z e= ^ 

We can, without loss of generality, choose this constant to be unity. Since u and there- 
fore x ^ nonconstant, the transformation 

z - X 

maps the region occupied by the body of water onto either the interior or the exterior of 
the unit circle of the x-pl^e- By replacing u by -u and hence w by -w, if neces- 
sary, we can further suppose with no loss of generality that the mapping 

z - X 

transforms the region occupied by the body of water onto the interior of the unit circle. 
Notice that ( (o is mapped onto the circle. The region occupied by the body of water in 
the x-pl 3116 is shown in figure 1. The existence of this function is a direct consequence 
of the Riemann mapping theorem which shows that every simply connected region bounded 
by a closed curve ^ can be mapped conformally into the interior of the unit circle. In 
addition, this mapping can always be chosen so that any given point in the interior of r € 
maps into the origin of the unit circle and any direction through the point corresponds to 
the direction of the positive real axis. Once this is done the mapping is unique. In fact, 
if 


z- Xl 

is any mapping of the interior of r £ into the interior of the unit circle, and if Zq is any 
point in the interior of c €, then the mapping z — defined by 
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X 2 (z) = 


X x (z) - Xx( z 0 ) 
1 - xJ(z 0 )Xi (z) 


also transforms the interior of ^ into the interior of the unit circle but in this case the 
point Zq maps into the center of the circle. 

We shall now show that once the boundary value problem for determining the surface 
displacement is transformed into the x~Pl ane it can be reduced to solving an ordinary 
differential equation. 

It has been shown that the region occupied by the body of water maps into the interior 
of the unit circle in the x~pl anc and that the shore line maps into the unit circle itself. 

We have also shown that with the change of dependent variable in equation (60) and with 
the new independent variables taken to be the polar coordinates (r and 8) in the x - ?!^ 6 , 
differential equation (50) for the surface displacement transforms into equation (64) or, 
when h/d > 1/2, into the approximate equation (65). In order to completely transform 
the boundary value problem for determining the surface displacement into the x^Pl^^ 
it remains only to transform boundary condition (40) into a boundary condition on the 
unit circle in the x - Pl ane - 

To this end notice that, since u is constant on $, the normal n to ^ is given by 


|Vu| 

Hence, it follows from the Cauchy-Riemann conditions that the components of n can be 
written as 


n — _ 

1 

O e (dw\* = . 

1 

n l " 

dw 

\dz / 

dw 


dz 


dz 


/dw^ dw dz 


and 


n 2 


dw 

dz 


-9*h 


\dz / dz dw 


dw 


dw 

dz 




dz 
) — 

dw 




( 66 ) 


/dw\* _ 

dw 

W/ 

dz 




Thus, it follows from equations (39) and (49) that the boundary condition (40) can be 
written as 
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where we have put 




{[■ 


a + H(E + iF) ( & + 
\3x 



= 0 


a e a* 1 ) ♦ to® 


( 67 ) 


Upon using the Cauchy- Riemann equations and the chain rule for partial derivatives this 
becomes 


Re 

k- 

-) + H(E + iF) /— + i — ) 


\4 

w/ \9u dv J 


or 


HF^i- 

3v 

When h/d > 1/2, figure 3 shows that 


HE ^ = R& 

9u 

this becomes approximately 



{ H-i)K + K =f } e 

dv 3u 

It now follows from equation (60) that 

HF^- HE^ + i®Ev = 

3v 9u 2 du 



V must satisfy the boundary condition 

V sinh H - sin H 
cos H + cosh H 



Finally, upon introducing the new independent variables r and 6 given by equation (22) 
and recalling that the boundary ^ maps into the unit circle in the x - pl ane > the boundary 
condition for V on the unit circle becomes 

pV r + qV^ + sV = y (6) for r = 1, -7r < 6 < n (68) 


where 
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p = sinh H(l) - sin H(l) 
cosh H(l) + cos H(l) 

q = H(l) - H(l) + sin H(l) 
cosh H(l) + cos H(l) 


sin H(l) sinh H(l) dH 
[cos H(l) + cosh H(l)] 2 dr 


r=lj 


r(9) = 




sinh H(l) - sin H(l) 
cos H(l) + cosh H(l) 



- 

/ _ \ 

* 

Re 

\dw/ 

r=l_ 




( 69 ) 


(70) 


When h/d > 1/2, equation (68) can be replaced by the approximate boundary condi- 
tion 


P? r + q? fl = y(0) for r = 1, -7T < 6 < n, - >- 
r 6 d 2 


(71) 


where now 


P = i ' 
q = H(l) - 1^ 


(72) 


and 


y(9) = R.& 

oil, &)(—)* 



\dw/ 

r=l 


(73) 


We therefore now show that the surface displacement for a sea or lake of arbitrary 
shape can be found by solving differential equation (64) in the interior of the unit circle 
subject to boundary condition (68). Or, if h/d > 1/2, the surface displacement is ap- 
proximately obtained by solving equation (65) in the unit circle subject to boundary con- 
dition (71). 
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SUMMARY OF EQUATIONS IN X" PLANE 


Upon comparing equations (64) and (68) (or when h/d >1/2 the approximate eqs. 
(65) and (71)) with equations (24) and (25), we see that this boundary value problem is a 
special case of the one solved in the preceding section. Hence, it follows from equa- 
tions (28) and (60) that the surface displacement is given by 


ICO 


cosh H + cos H V* 


g = “ ~ r (r)e 

sinh H - sin H /—d n 

n=-°° 


in0 


or, when h/d > 1/2, approximately by 


5 - fi " (r)el 

n=-oo 


1110 for *>± 
d 2 


(74) 


(75) 


where (as can be seen from eqs. (22), (32) to (34), (56), and (57)) each function 
O n (r) is determined by solving the ordinary differential equation 


r(rS2 


<■>' ' (’ 


r — + a; 2 + n 2 - 2in/3)f2 
dr 


r_ < 


n = 0, ±1, ±2 . . . , 
0 < r < 1 


( 76 ) 


subject to the boundary conditions 


P«„(l) + (s + inq)fi n (l) = y 


n' 


lim S2 (r) is finite 
n 


n 


> n = 0, ±1, ±2, . . . 


r— 0 


( 77 ) 


The right sides of equations (76) and (77) are given by (see eqs. (31), (63), and (70)) 



r 1 /cosh H + c os H f e _ in0 H 

n 2? r V sinh H - sin H / L 

*'-77 




[>> + 

4 2) 1 

dz 

L x 

y j 

dw 


d<9 


1 /si 
7n 277 Vco 


sinh H(l) - sin H(l) 
cosh H(l) + cos H(l) 


I. 


f n = 0, ±1, ± 2 , . . 


-ine 


Re 


r / v* 1 

oil, e) f 

\dw/ r= i 


ae 


(78) 


or, if h/d > 1/2 (see eqs. (31), (65), and (73)), these can be replaced by the approxi- 
mate relations 



The coefficients a and 0 are given by (see eqs. (51) and (22)) 


(79) 


sin H sinh H „ dH 

a = - — r — 

(cosh H + cos H)(sinh H - sin H) dr 
. 1 sinh H + sin H dH 

/3 = r — 

2 cosh H + cos H dr 


(80) 


and when h/d >1/2 these can be replaced by the approximate relations (see eqs. (22) 
and (53)) 


a ~ 0 


*\ 


2 dr 


). for 


(81) 


Finally, the constants p, q, and s are given by (see eq. (69)) 
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_ sin h H(l) - sin H(l) 
cosh H(l) + cos H(l) 

q = H(l) - s inh Hi 1 ) + sln H ( 1 ) 
cosh H(l) + cos H(l) 


"\ 


( 82 ) 


s = - 


sin H(l) sinh H(l) dH 
[cos H(l) + cosh H(l)] 2 dr 


r=l 


or when h/d > 1/2 approximately by (see eq. (72)) 


P = 1 

q = H(l) - 1 l 
s = 0 


(83) 


Thus, an exact solution for the surface displacement (and since all other physical 
quantities can be expressed in terms of this, an exact solution to the complete flow prob- 
lem) can be obtained once the ordinary differential equation (76) with boundary condition 
(77) has been solved. This solution is given parametrically in terms of the variables r 
and 6. These parametric variables are determined in terms of the physical variables 
x and y by the mapping 


z — X 
iG 

with z = x + iy and x = re > which maps the lake or sea conformally into the unit 
circle with a selected point being mapped into the origin. The coefficients of the ordi- 
nary differential equation depend upon the particular choice of the function H(r). Thus, 
it is impossible to proceed further in the general case. 

The procedures involved in obtaining a complete solution are best illustrated by 
considering a particular case. In order to simplify the calculations it will be assumed 
that h/d >1/2 and approximate forms in equation (81) will be used for the coefficients 
a and /3. A particular depth distribution function H of r will first be chosen. The 
ordinary differential equation will then be solved and the surface displacement will be 
found as a function of r and 6. A particular shape for the body of water will then be 
chosen. Once this is done the relation between the physical variables x and y and the 
parametric variables r and 0 will be found. This will give the surface displacement 
as a function of x and y. 
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SOLUTION OF ORDINARY DIFFERENTIAL EQUATION FOR 
PARTICULAR BOTTOM TOPOGRAPHY 


We shall consider the case where the function H of r is given by 


H(r) = H + 6(1 - r 2 ) 0 < r < 1 

o — — 


(84) 


where H„ is the constant depth at the shore line and H + 6 is the depth of the deepest 
s s 

point of the lake. Upon inserting equation (84) into equations (81) and (83) and then in- 
serting the results into equations (76 and (77), we obtain the following boundary value 
problem for 


r 2 n'' + rS2’ - fn 2 + 2in6r 2 ~h = T 
n n |_ J n n 


n = 0, ±1, ±2, . . . 
0 < r < 1 


(85) 


n^(i) + m[H s - i]syi)=r n 

lim ^ n (r) is finite 


n = 0, ±1, ±2, . . . 


( 86 ) 


Before obtaining the solution to this problem, we shall first show that 



~ r o( r)dr = y 0 

r 


(87) 


To this end notice that equations (79) imply 



J 


Hence, it follows from the Cauchy-Riemann equations and equations (18) and (22) that 
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27 r 


277 


277 



4 2) 1 

/dx dy 

- de 

y 

\0u dv 

dv 0u/r 

ofl 

/dx dy . 

. 9x9y\ drd0 


\0r 00 

00 0r/ 

4 2 >! 

a (*>y) , 

dr d0 

y J 

0(r, e ) 



Upon recalling that the unit circle in the x - Pl ane transforms into the interior of the body 
of water in the z-plane, it follows from the theory of transformation of multiple inte- 
grals that 



r Q (r)dr 




dx dy 


(89) 


where the integration is carried out over the entire surface of the region occupied by the 
body of water. Let S denote the distance measured along the boundary ^ of this re- 
gion. Then applying the Divergence Theorem to equation (89) yields 



- T 0 (r)dr 
r 



+ CT (2) n 2 dS 


Hence, upon substituting equations (66) into this relation, we obtain 
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Now 


and 


Hence, 




dS =\l( dx) 2 + (dy) 2 


dw 

dz 


& j(f ) - ’fa* - <»/ 


dw 


dz 


ds = y[ u x dx j + [ v x dy ] 2 + [jx dx ] 2 + 1 x dy j 


Upon using the Cauchy- Riemann equations (11) and noting that u u = -vv , we find that 

x y y x 


dw 


dz 


dS = t/ (u v dx + u y dy) 2 + (v x dx + v y dy) 2 = J{ du) 2 + (dv) 2 




Hence, it follows from equations (22) 



But r is constant (and equal to unity) on Hence, 


dw 

dz 


dS = d6 


and equation (90) becomes 



- T 0 (r)dr 
r 


27 T 



Re 


/dz \* 

o(l, 9) I — ) d 9 

\ dw / r= l_ 


Finally, comparing this with equation (88) shows that equation (87) holds. 

Having established this result, we can now proceed to solve boundary value problem 
determined by equations (85) and (86) for n = 0. In this case equations (85) and (86) become 

34 


■mil 


■ min i m mi i ii i 


I llllll INI 


l 


r(rOj,)' = r 0 0 < r < 1 


^(1) = y 0 


lim ^(r) is finite 


Upon integrating equation (91) we get 


r^Q = Zj + 


f * 0 

y 0 


P 


where is a constant of integration. Upon substituting this expression into the first 
boundary condition (eq. (92)), we find, upon using equation (87), that l j = 0. Therefore, 


% • r r o 

*4) 


<p)^ 

p 


Integrating this with respect to r yields 


2 o o f | f - r 0 


(p)dp dr 


where Iq is a constant of integration. Upon integrating by parts, we find that 


%(r) = l Q + In 


r J 

r f r 0 (p) + f in p r«(p) 

Jo p J v p 


Hence, upon defining the function of r and p by 


<#0 (r » & = 


i In r for r > p 


In p for p > r 
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the solution ^(r) can be written as 


sy r ) - + 71 f M 0 <?iP) T 0<P) (94) 

Jq p 

Notice that this solution automatically satisfies the second boundary condition (92). 
This would not be the case if equation (87) did not hold; for then l ^ would not vanish 
and the solution would therefore contain a term l ^ In r which is unbounded at the 
origin. On the other hand, the solution given by equation (94) is indeterminate to within 
an arbitrary constant Iq. But this is as it should be since the n = 0 term in equation 
(75) is independent of 6, and, as can be seen from equations (65) and (71), the original 
boundary value problem for £ involves only the partial derivatives of £ and hence 
determines £ only to within an arbitrary constant. 

We now return to the boundary value problems determined by equations (85) and (86) 
and suppose that n / 0. Put 


= y-2in 5 n = ±1, ±2, . . . 

r n s V n = **■> • • • 


(95) 


Then equation (85) becomes 


2 

a d £2 / 9 a\ 

2 + r n — - l” 2 - r n) n n = F n n = ±1, ±2, 

n 2 n d V n/ n n 
dr n 


dft. 

dr 


But this is just the nonhomogeneous Bessel’s equation which has the solution 


= O" + 
n n n 


(96) 


where and are the homogeneous and particular solutions, respectively. The 
homogeneous solution is 


- i i: n) w + 4 n) w ” ■= * 2 — 


(97) 


where for i = 1, 2 are arbitrary constants and and Y n are the n^* 1 - order 
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Bessel functions of the first and second kind, respectively. The Wronskian of J n (r n ) 
and Y n (r n ), denoted by W(J n , Y n ), is 


W(J„, Y„) . 


J (r ) Y (r ) 
n v n 7 n v n 


dJ„( r _) dY_(r ) 
n' n 7 n v n 7 


div dr 
n n 


7rr 


n 


Hence, a particular solution £2^ of nonhomogeneous equation (85) is (ref. 7, p. 155) 


n £ = - 
n 2 


[* J n<V> Y n<V> r n<P> ^ + / WW r nW 
Jn P Jr P 


n = ±1, ±2, . . . 

(98) 


This solution is bounded. To show this, notice that (ref. 9, pp. 15 and 62) for r — 0 

"\ 


n 


J n<V> 


I ±— 

on, 


H- 


r n = ±1, ±2, . . . 


J 


Hence, by introducing these asymptotic forms for and Y n into equation (98), we find 
that 4 


lim n£(r) = — T n (0) n = ±1, ±2, . . . 
r—0 [ n | 2 

Thus, fi^(r) and J n (X n r) are bounded at r = 0 and Y n (X n r) is unbounded. Hence, equa- 


4 It is not hard to verify that r n (0) = 0 for all n. Therefore, 

lim fl^(r) = 0 for n = ±1, ±2, . . . . 
n ’ ’ 

r— 0 
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tions (96) and (97) show that the solution S2 n is bounded if, and only if, = 0. The 
second boundary condition (86) can therefore only be satisfied if 


Z< n) = 0 for n = ±1, ±2, . . . 
Upon combining equations (96) through (99) we get 


(99) 


n n< r) = i l” )j n< X n r > + ! 


r J n<V> Y n< x n r > r n<P> ~ * /' J n< x n r > W> r n<<» 
Jf\ P Jf 


dp 


( 100 ) 


Hence 


!J n (l) = 4 n) J„(x n ) + | Y n (X n ) f 1 J n (x n p)r n 
• 4 n)j n< x n> + f Y X> f J n<V> r n 


d£ 


<P>^ 

P 


where J^(X n ) = j^(d/dr)J(X n rjJ ^ and so forth. Substituting these results into the first 

boundary condition (86) and using the recurrence relations for the Bessel functions to 
eliminate their derivatives gives 


r (n) 1_ 


CD, 


( 1 ) 

n *- 


^ y n p 


(101) 


where 


"n ’ = n [f + ‘< H S - 

‘ " [l + ‘< H s - Y n<V " X n Y „ + l< X n> 
Inserting equation (101) into equation (100) shows that 


> n = ±1, ±2, 


(102) 
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fi n (r) = -^77 J n<V> + f f ^ r n^ ~ n = 

o> (1) 2 •/() P 


= ±1, ±2, 


where we have put 


cu 


(2) 


Y ” (x ” r) - Jr, w> 

n 




J n (x n p) 


for p < r 


^ n (r; p) =-{ 


V n = ±1, ±2, 


OJ 


(2) 


Y „<V> - J n<V> 

C X ) ' ' 

n 


J „<V> 


for r < p 


Since 




and 


Y -„< X -„ r > = (-‘'"[VV*]' 

and since equations (79) show that 


r = r* 

-n n 


y-n^n 


it is easy to see from equations (102) to (104) that 


p -l¥ 

o Pl = (-l) n 4 J) for j = l > 2 


and hence that 


n (r; p) = jj# n (r; pfj 


(103) 


(104) 
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Therefore. 


fi_ n = for n = ±1, ±2, . . . (105) 

Upon substituting equations (103) and (94) into equation (75) and using equation (105), 
we find that the solution for the surface displacement is 


?(r, 9) = l Q + 2/?e j J n ( x n r ) eUiy + * 


in (9 





? in0 f l ,y{ n (v. p )T n {p) 

Jq 


dp 

p 


n=l n-0 (106) 

where, upon collecting the definitions of the various quantities in this equation, we get 


r i 


— In r for p < r 
7 r 


M 0 (r; p) = J 




n- 


•^ n (r; p) =•< 


u>. 


( 2 ) 


Y (X r) - — ?L_ J (x r) 
n v n 7 n v n 7 

n 


4 2) 

W> ' W> 

ur 7 

n 


— In p for r < p 

7T 


J n (\p) for p < r 


J n (*„r) for r < P 


J 


(93) 


> n = 1, 2, . . . 


(104) 


«?> - "[ 

} + ‘< H S - ‘OW - Vn +l <V 


i + i(Hg - D]Y n <V - X n Y n+1 (X n ) 


X n = y~2in 5 = e 377i / 4 ^2n6 = (i - 


> n = 1, 2, . . . 


(102) 


(95) 


40 


ia I I 


r »^ / Vta 1 

•'-IT 


foM + </ 2) l 

dz 

L x y J 

dw 


d9 


> n = 0, 1, 2, 


= JL e -ine 

y ° 2 tt / 


Re 


a(l, 0) 

\dw. 


r=l 


d0 


o.^ + to* 2 ) 

Finally, upon inserting equation (49) into equation (47), we find that 


<,<» = rJe- H /2 2 sin ~ + rj (l - e" 0 / 2 2 cos ^ 
a (2) = T t^-H/2 2 cos H.i) +T t e -H/2 2sin H 


in which 


and where we have put 


H(r) = H g + 6(1 - r^) 0 < r < 1 


t 2lTT i 

r. = for i = 1, 2 

1 gd 


(79) 


(67) 


(107) 


(84) 


(108) 


If h/d is larger than 1, then equations (107) may be replaced by the less exact re- 
lations 




(48) 


Notice that when these latter approximate relations are used for cr^ and and 
when the components of the components of the wind stress and Tg are constants, then 
(see eq. (79)) 

r n = 0 for n = 0, 1, 2, . . . 
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and only the first summation appears in solution (106). 

In order to evaluate the integrals in equation (79), it is necessary to know the rela- 
tion between the physical coordinates x and y and the polar coordinates r and 9 in 
the x~pl ane > and also to know dz/dw = x(dz/d\) as a function of x (or equivalently as a 
function of r and 9 ) . These relations can be determined as soon as the function 

z - X 

which maps the interior of the region occupied by the body of water in the physical plane 
into the interior of the unit circle in the x - plane is selected. 


SOLUTION OF COMPLETE BOUNDARY VALUE PROBLEM FOR ELLIPTICAL LAKE 

OR SEA BY CONFORMAL MAPPING 

In order to find the mapping z — x, the specific shape of the body of water must be 
chosen. We shall, therefore, suppose that the body of water has the shape of an ellipse 
whose semimajor axis has length and whose semiminor axis has length Lg as 
shown in figure 6. In addition, we shall assume that the body of water is deepest at its 
center. Since the depth-distribution function (eq. (84)) has its maximum at the center of 


y 
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the unit circle, the mapping z — x must be chosen to map the center of the ellipse into 
the center of the unit circle. For convenience the mapping will be chosen so that it maps 
the positive real axis in the z-plane into the positive real axis in the x~Pl ane - As has 
been shown in the preceding discussion, these conditions are sufficient to completely 
determine the mapping 


z — X 

In order to introduce this mapping it is convenient to introduce the parametric com- 
plex variable 


t = ? + ir? 


(109) 


The t-plane is shown in figure 7 . Let 


K = K(k) 

be the complete elliptic integral of the first kind of modulus k and let 

K’ = K(k') 


where 



- k 


2 


( 110 ) 


(HI) 


( 112 ) 
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Figure 7. - t-plane. 
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is the complementary modulus. Then the mapping 


z = I sin — 
2K 


(113) 


where 


I s 




(114) 


transforms the rectangular region in the t-plane shown in figure 7 into the elliptical re- 


y 



gion (with two slits) in the z-plane shown in figure 8. The mapping 

X=k*^sn(t, k) (115) 

transforms the rectangular region in the t-plane shown in figure 7 into the interior of 
the unit circle with two slits in the x _ Pl ane shown in figure 9. Hence, the combination 
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Figure 9. - X'Plane. 

of mapping (115) with the inverse of mapping (113) is the mapping (ref. 10, p. 371) 



(116) 


which transforms the elliptical region with two slits in the z-plane (shown in fig. 8) into 
the interior of the unit circle with two slits in the x - Pl ane (shown in fig. 9). Since each 
cut in the z-plane corresponds precisely to a cut in the x~Pl ane , these cuts can be sealed 
and hence, mapping (116) is the one with the desired properties. However, it is more 
convenient to deal with this mapping in terms of its components (mappings (113) and 
(115)). It is easy to verify that the modulus k is determined in terms of the physical 
dimensions of the ellipse by Jacobi's nome (ref. 10, p. 371) 


e-iK'/K 



(117) 


Upon taking the real and imaginary parts of equation (113), we find that the physical 
variables x and y are given in terms of the parametric variables ij and r\ by 
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= I sin cosh r}j 


y = I cos [ — Asinh (— ri 
\2K ) \2K j 


Upon taking the real and imaginary parts of equation (115) and recalling that 


(118) 


re 10 = x 


we find that the polar coordinates r and 9 in the x _ P lane can be expressed in terms of 
the parametric variables ij and t] by (ref. 11, p. 24) 


tan 6 = 

cn(£, k)dn(£, k) 

sn(? 7 , k’)cn(? 7 , k’) 


sn(£, k) 

dn(? 7 , k') 


r cos 6 = 


y k sn(g,k)dn (77, k') 

_1 - sn 2 (? 7 , k’)dn 2 (5, k)_ 


(119) 


These equations now determine the physical coordinates x and y parametrically (with 
parameters £ and r\) in terms of the polar coordinates r and 6 in the x _ Pl ane - Actu- 
ally it is possible to eliminate the parametric variables | and 77 between equations 
( 118 ) and (119) and obtain an expression for x and y directly in terms of r and 6 but 
we shall not do this here. 

Finally, in order to evaluate the integrals in equations (79), it is necessary to obtain 
an expression for dz/dw. Because equation (20) implies 


dz _ ^dz 
dw dx 


it follows that 


dz _ ^dz d^ 
dw dt dx 

Upon differentiating equations (113) and (115), we find that 
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and 


Therefore, 


dt = 1 

dx k 1//2 cn(t, k)dn(t, k) 


dz 

dw 


fflX 

2^k K 



cn(t, k)dn(t, k) 


(120) 


In order to separate the real and imaginary parts of the denominator of this expression, 
notice that the addition theorems for elliptic functions (ref. 11, p. 23) show that 


sn(t + t*) - sn(t - t*) = 2 sn t* cnt _dnt_ 

9 9 9 

1 - k^ sn t sn t* 


Hence, because sn(t*) = (sn t)*, we see that 


sn(2(j) - sn(2irj) = 


2(sn t)* cn t dn t 


1 - k 


2 


sn t 


4 


It now follows from equation (115) that 

1 = _2 X* 

cn t dn t ^k (l - | x | 4 j |sn(2£,k) - sn(2irj, k)J 

= X* 1 

^k (l - | x | 4 ) sn(2 ^ k) ’ itn (27?,k') 

(In more recent references the symbol sc is used instead of the symbol tn used 
here.) Substituting this into equation (120) and using equation (21) yields 

„ cos (— t) 

jlz _ jrf \2K / 

dw kK j _ r 4 sn(2£,k) - i tn(2rj, k') 
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or 


/ 2 \ cos (— |Vosh/ r — rj\ - isinf— lAsinh f—r]\ 

dz _ 7r I | r 2 \ \2K / \2K ) Uk / \2K ) 

dw kK _ r 4 y sn(2|, k) - i tn(277, k’) 


0 <r < 1, -7T < 0< 3T 


( 121 ) 


This expression is indeterminate on the boundary of the unit circle (r = 1). It is, there- 
fore, necessary to obtain an alternate form for this expression on the boundary. To this 
end notice that figure 7 shows, at r = 1, r\ = ±(K*/2). 

Now it follows from the addition theorems for the elliptic functions that 


cn 


'{ * 1 £.)dn ± i iiis c n | dn j - k sn 2 { 

> 2 / l 2 / y k 


cn £ dn £ t i(sn £ dn 2 i; + k cn 2 | sn |) 
(1 + k sn 2 £) 2 


Since (ref. 11, p. 20) 


dn 2 £ = 1 - k 2 sn^ £ 


and 

9 9 

cn | = 1 - sn £ 


this becomes 


, K 


cn [I ± i — )dn ( £, ± i 


, K* 


1 + k / l - k sn 2 | 
^k \1 + k sn 2 £ 


cn i dn * t (i + k)i — ^nj — 
1 + k sn 2 £ 1 + k sn 2 £ 


However, equations (119) implies that at r = 1 (or rj = ±(1/2)K*) 

cos e = I* + k)sn 1 " 

1 + k sn 2 £ 

► 

sin g = ± cn | dn g 
1 + k sn 2 | 


( 122 ) 
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Hence, 


cn 


(* 



Upon inserting this into equation (120) evaluated at r = 1, we find that 


dz 

dw 


= ±1 


jrl /l+ksn ^ 1 

r=l 2(1 + k)K L . k sn 2 


Icos 


' I ±i» K ° 


t2K 


4K 


_ itI /l + k sn^ ; 

‘ 2(1 * k)K ^ . k sn 2 { 

5 

But it follows from equations (117) 
bolic functions that 


\ sin (JL A sinh (* ± i cos (JL A cosh 

/ \2K / \4 K/ \2K / \4K/ 

and the use of the exponential forms of the hyper- 


and 



5 Note that 


and 




- L 


2 


+ L 


2 
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Hence, 


(123) 


dz 

dw 


r=l 


1 + k sir g 
2 < 1 + W K Vl-ksn 2 {, 




where the plus sign corresponds to 9 > 0 and the minus sign is used when 9 < 0. Since 
equation (123) is to be used in the integral of the second equation (79) and since the inte- 
gration here is with respect to 9, it is convenient to have an expression for the variable 
5 as a function of 9. This can be accomplished by solving equation (122) for sn (j as a 
function of 9 to obtain 


sn £ 


1 + k ±yi + k 2 - 2K cos 29 
2k cos 9 


Since £ = K at 9 = 0, it is easy to see that the minus sign must hold in this equation 
and therefore that 


£ = sn 


-1 


1 + k 




k - 2k cos 29 


2 k cos 9 


,k 


for r = 1 


(124) 


Since 77 = ±(K'/2) when r = 1, it follows from equation (118) and the relations imme- 
diately preceding equation (123) that 


x = L, sinf— % 
L \2K 


r for r = 1 




(125) 


Thus when the wind stress components and are given as functions of x and y, the 
integral in the second equation (79) can easily be evaluated by using equations (67), (84), 
(107), (108), and (123) to (125). Similarly the integral in the first equation (79) can be 
evaluated by using equations (84), (107), (108), (118), (119), and (121). These results 
can then be substituted into equation (106) to determine the surface displacement £ as a 
function of r and 9. Equations (118) and (119) will then give £ as a function of x and 
y and this completes the solution to the problem. 

Since the geometry chosen does not represent any particular lake or sea, it will not 
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be worthwhile to carry out a complete parametric study for various functional forms of 
the wind stress. Instead we shall only obtain numerical results for the simplest case 
which illustrates all the features involved. Thus, we shall suppose that h/d > 1 and 
that approximate equation (49) can be used for and o^\ We shall also sup- 
pose that the wind stress is constant and is given in dimensionless form by 

t t 

T i = Tq cos 3 



t 

T 0 


sin 3 


> 


(126) 


where Tq is the constant dimensionless magnitude of the wind-stress and 3 is the angle 
between the wind stress and the x-axis. Then in this case, as indicated previously, 

T n = 0 for n = 0, 1, 2, . . . and it follows from equations (48), (67), (79), (106), (108), 
(123), and (124) that the solution for the surface displacement is 


s - ? 


o 


t L 2(1 + k)K 
0 1 


Re 


E — U* 

n=l 


(127) 


where 



DISCUSSION 

A method for obtaining solutions to a certain class of partial differential equations is 
introduced. This class of equation includes Welander's equation for the wind-driven cir- 
culation in shallow seas and lakes for a large class of bottom topographies. The tech- 
nique, together with some results based on conformal mapping is used to reduce the prob- 
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lem of finding solutions to Welander's equation for a closed body of water to the problem 
of solving an ordinary differential equation which depends on a parameter . An approxi- 
mate form of this ordinary differential equation which is valid when the depth of the body 
of water is larger than one half of the Eckman depth is solved for a particular bottom 
topography. This solution takes on a particularly simple form when the wind stress is 
constant and the depth of the lake is greater than the Eckman depth. 

The following discussion concerns the application of this solution to an elliptically 
shaped body of water with an approximately paraboloidal bottom topography whose depth 
is larger than the Eckman depth and for which the wind stress is constant. No attempt 
is made to conduct a parametric study to demonstrate the effect of the governing dimen- 
sionless parameters since this would be of practical value only if a particular lake or sea 
were being studied. For this reason numerical results will only be given for the surface 
displacement of the surface of the body of water. All other flow quantities of interest can 
easily be expressed in terms of this quantity by using the results given in reference 1. 

The surface displacement of the body of water is given in terms of the intermediate 
variables r and 9 by equations (127) and (128) where the modulus k of the complete 
elliptic integrals is determined from equation (117) and the parameters u>^ and X n are 
given by equations (102) and (95), respectively. Since the physical coordinates x and y 
of the horizontal plane of the body of water are related to the intermediate variables r 
and 9 parametrically by equations (118) and (119), the surface displacement is known as 
a function of x and y. These results involve the physical parameters Lg/Lj, H g , 6, 
and 3. The first three of these parameters characterize the geometry of the body of 
water and the last one is the direction of the wind stress. All the numerical calculations 
were performed for H = 5 and 6 = 4. The results are shown as isometric projections. 
The results of the calculations for 3 = 0 (i. e. , when the wind direction is along the ma- 
jor axis of the ellipse) are shown in figure 10. The figure shows that the water piles up 
downwind and the surface is s-shaped across the wind direction. The results for 
3 = 45° and 3 = 90° are shown in figures 11 and 12, respectively. In all three figures 
the ratio Lg/Lj of the semiminor to the semimajor axis of the ellipse is taken as 2/3. 
The effect of fixing 3 at 45° and changing L^/L^ to 1/2 is shown in figure 13. It can 
be seen from these figures that the amount of surface displacement is essentially deter- 
mined by the distance across the lake in the direction of the wind. The greater this dis- 
tance, the larger the surface displacement. 
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Figure 12. - Surface displacement of an elliptic lake. Ratio of minor to major axis, LJL, ■ 2/3; wind parallel to minor axis 
(o = 90°); 8=4; H s = 5. 



Figure 13. - Surface displacement of an elliptic lake. Ratio of minor to major axis, LJli ■ 1/2; wind at 45° to axes 
(J =45°); 6=4; H s = 5. 


CONCLUSIONS 


A technique for obtaining solutions to a certain class of partial differential equations 
is given. This technique is applied to the equation which describes the wind- driven cir- 
culation in a shallow sea or lake, ft is shown that for a large class of depth distributions 
this technique can be used to reduce this partial differential equation to an ordinary dif- 
ferential equation. Complete analytical solutions are obtained for a specific geometry of 
the lake or sea. 
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